Turbulent transport in tokamak plasmas with rotational shear 
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Nonlinear gyrokinetic simulations have been conducted to investigate turbulent transport in toka- 
mak plasmas with rotational shear. At sufficiently large flow shears, linear instabilities are sup- 
pressed, but transiently growing modes drive subcritical turbulence whose amplitude increases with 
flow shear. This leads to a local minimum in the heat flux, indicating an optimal E x B shear value 
for plasma confinement. Local maxima in the momentum fiuxes are also observed, allowing for the 
possibility of bifurcations in the E x B shear. The sensitive dependence of heat flux on temperature 
gradient is relaxed for large fiow shear values, with the critical temperature gradient increasing 
at lower flow shear values. The turbulent Prandtl number is found to be largely independent of 
temperature and flow gradients, with a value close to unity. 



Introduction. Experimental measurements in mag- 
netic confinement fusion devices indicate that sheared 
mean E x B flows can significantly reduce and some- 
times fully suppress turbulent particle, momentum, and 
heat fluxes [1] [2] . Since these turbulent fluxes determine 
mean plasma density and temperature profiles, their re- 
duction leads to a local increase in the profile gradients. 
This increase can be dramatic: transport barriers in both 
the plasma core and edge have been measured with radial 
extents on the order of only tens of ion Larmor radii [3] . 
The associated increase in core density and temperature 
results in increased fusion power. Thus, understanding 
how shear fiow layers develop and what effect they have 
on turbulent fiuxes is both physically interesting and 
practically useful. 

This Letter reports a numerical study of the infiuence 
of sheared toroidal rotation on turbulent heat and mo- 
mentum transport in tokamak plasmas. Two main ef- 
fects of sheared toroidal rotation were identified in previ- 
ous numerical work W^ : suppression of turbulent trans- 
port by shear in the perpendicular (to the mean magnetic 
field) velocity and linear destabilization due to the par- 
allel velocity gradient (PVG). While the former observa- 
tion indicates that a finite fiow shear improves plasma 
confinement, the latter raises the question of whether 
more shear is always beneficial. Below we report that 
the PVG-drivcn linear instability |9j is stabilized at suffi- 
ciently large flow shear values, consistent with fluid the- 
ory in slab geometry [10]. Correspondingly, fluxes de- 
crease with increasing flow shear as the linear stabiliza- 
tion point is approached. However, beyond this point, 
transiently growing modes driven by the PVG give rise 
to subcritical turbulence. The fluxes associated with this 
turbulence increase with flow shear. This implies an op- 
timal flow shear for each temperature gradient; the fact 
that the minimum heat flux value is finite indicates that 
there is a maximum attainable temperature gradient that 
can be maintained for a given heat flux. Additionally, the 
observed presence of maxima in the momentum fluxes ad- 
mits the possibility of bifurcations in the flow shear (and 
thus the temperature gradient). 



In the absence of flow shear, a small increase in temper- 
ature gradient leads to a large increase in heat flux ( "stiff 
transport"). Recent experimental evidence |11) suggests 
that flow shear may reduce this sensitivity in configura- 
tions with low magnetic shear. Our results indicate that 
at low flow shear values, both the critical temperature 
gradient for the onset of turbulence and the stiffness in- 
crease. At high flow shear values, the opposite behavior is 
observed (the stiffness and critical temperature gradient 
both decrease). 

Model. A closed set of evolution equations [121 [13] 
for mean plasma density, pressure, and toroidal angu- 
lar momentum is obtained by taking moments of the 
kinetic equation and applying the 5f gyrokinetic order- 
ing |14Lll5j . These moment equations relate the evolution 
of the mean quantities to particle, momentum, and heat 
fluxes, which are typically dominated by turbulent con- 
tributions. We restrict our attention to electrostatic fluc- 
tuations and assume a modifled Boltzmann response |16| 
for the electron distribution. The resultant particle flux 
is identically zero, and the radial components of the tur- 
bulent heat flux of species s, Qa, and toroidal angular 
momentum flux, 11, are given by 
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where ■0 is a flux-surface label, (j) the toroidal angle, rus 
the particle mass, v its velocity in the frame of mean 
flow, R its major radius, v^; the fluctuating E x B drift 
velocity, Sfs the deviation of the distribution function 
from a local Maxwellian, and the overline denotes a spa- 
tial average over a thin annular region encompassing a 
given flux surface. 

The distribution function Sfs appearing in Eqs. M 
and (pi) is calculated by solving the standard Sf gy- 
rokinetic equation in the limit where the plasma flow 
speed, u, is ordered comparable to the ion thermal speed. 
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FIG. 1: Average linear growth rates vs. flow shear for differ- 
ent values of Rq/Lt- Inset: critical flow shear vs. Rq/Lt- 



■i^th [131 HZ] ■ In this "high- flow" regime, the flow veloc- 
ity is constrained to be u = R'^uj{tp)\/(j), where R is the 
major radius and uj is the rotational frequency |18| . We 
consider a subsidiary expansion in low Mach number M 
{pi/L -^ M <^ 1), where L is a typical macroscopic 
length scale (e.g., the minor radius of the torus) and pi 
is the ion Larmor radius. However, we allow for flow gra- 
dients of order 1/ML, and so neglect terms proportional 
to Lo (such as Coriolis and centrifugal drifts) but retain 
those proportional to du/d-ip (flow shear). The resulting 
equation is 
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where d/dt = d/dt + u-W, h = 6f + {e<f/T)Fo is the 
deviation of 6f from a Boltzmann response, e is particle 
charge, (p is the electrostatic potential, v^ — (c/i3^)B x 
V(^, T is temperature, Fq is a Maxwellian distribution 
of peculiar velocities, C is the collision operator, B is 
the magnetic field strength (with B^ the toroidal com- 
ponent), vs = h/il X (wyb • Vb -I- v']^S/B/2B] contains 

magnetic drifts, fi is the Larmor frequency, b = B/B, 
and (.)jj denotes the average over gyro-angle at fixed 
guiding center position R. 

Eq. p| is solved in the rotating reference frame using 
the local, nonlinear gyrokinetic code GS2 [8l [19]. The 
toroidal velocity is expanded in ip about its value at the 
center of the simulation domain, ■00 j giving u « R^{ip — 
'ipo){doj/d'tlj)V(j). The local approximation assumes that 
duj/dip is constant across the simulation domain. Con- 
sequently, the effect of sheared flow is accounted for 
by a single parameter, je — {i'/<l)[dijj/d'4')RQ/vth = 
{M /q){d\nuj/d\x\r), where q is the safety factor, Rq is 
the major radius at the center of the flux surface, and r 



is the half-diameter of the flux surface (both measured 
at the height of the magnetic axis). 

We study a system whose magnetic geometry corre- 
sponds to the widely-used Cyclone base case [20] (un- 
shifted, circular flux surface with q = 1.4, magnetic shear 
s = dhvq/dXnr = 0.8, r/i?o = 0.18, and i?o/L„ = 2.2, 
with L~^ = — din n/dr). The E x B shearing rate, 7^;, 
and the normalized inverse temperature gradient scale 
length, K, = Rq/Lt, were varied over a wide range of 
values in a series of linear and nonlinear simulations. 

Linear stability. The average linear growth rates, 7, 
obtained from these simulations are given in Fig. [l] We 
see that 7 decreases rapidly with 7^ for small values of je 
before increasing to a local maximum and subsequently 
decreasing to zero. For je ^ 0.25, the system is linearly 
unstable only in the presence of both temperature and 
parallel flow gradients (cf. [SHE]). The system becomes 
linearly stable for large values of je, with the critical je 

(c) 

for stability, 7^ , increasing approximately linearly with 

K. Beyond 7^ , there are no linearly unstable modes. 
This result is qualitatively similar to the prediction from 
fluid theory in a slab [10], where no linear instability is 
possible when 7b/s exceeds a certain critical value. 

Heat flux. The turbulent heat flux calculated from 
nonlinear simulations is given in Fig. [2] For all k val- 
ues, the heat flux qualitatively follows the same trend as 

(c) (c) 

growth rates when je < 7b • For 7^ > 7^; , nonlinear 
simulations initialized with low-amplitude noise develop 
no turbulent transport. However, for finite initial fluctu- 
ation amplitudes, one finds that the turbulence does not 
necessarily decay away: for sufficiently large values of k 
and initial amplitude, the flux reaches steady state values 
in excess of those found for je just below 7^^ ^T\. The 
flux then increases monotonically with je ■ For the range 
of jE considered here, the heat flux associated with a 
given temperature gradient is thus minimized at a finite 
shearing rate. 

Subcritical turbulence. Because the tiu'bulence is 
present in the absence of linear instability, we refer to 

(c) . . . 

it as subcritical. When 7^ > 7^ , linear simulations 
exhibit transient growth, with order unity increases in 
the initial fluctuation amplitudes over times of several 
-Ro/wth before subsequent decay (Fig. [3|. The duration 
of growth decreases with increasing flow shear, but the 
transient growth rate increases so that the amplification 
factor of the initial perturbation amplitude grows with 
flow shear (cf. [10 ). This transient amplification pro- 
vides an energy source for the turbulence, which can be 
maintained by the nonlinearity through redistribution of 
energy amongst other modes. We have identified the 
PVG term in Eq. [s] (the duj/dip term on the RHS) as 
the driver of the subcritical turbulence, as the latter is 
no longer present when the PVG is artificially set to zero. 
Momentum flux. The toroidal angular momentum 
flux also mimics the linear growth rate, except near 7^; = 
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FIG. 2: Turbulent heat (top) and toroidal angular momentum 
(bottom) fluxes vs. flow shear for various values of Rq/Lt- 
For 7_B < 0.25, the heat flux is reduced hy E x B shear. 
Beyond this value, the PVG drives instability and both heat 
and momentum fluxes increase. For 7_b > 7^ , the turbulence 
is subcritical, sustained by transiently growing modes. 







6 8 10 

t{vth/Ro) 

FIG. 3: Heat flux versus time obtained from linear simulations 
for Rf)/LT = 8.75 and three 7_b values at which subcritical 
turbulence is observed. 



FIG. 4: Turbulent Prandtl number as a function of flow shear 
for various values of Rq/Lt- Unfilled points correspond to 
linearly stable flow shear values. Zero Pr points represent 
fully suppressed turbulence. 



[22j . This is understood by expressing the momentum 
flux in diffusive form: 11 = —raiVt\i{qRo/r)vi{K,^E)lE, 
where Vi is the turbulent viscosity. For 7^; small, the 
fluctuation amphtudes vary little, so Vi is approximately 
constant and 11 oc 7^;. For larger 7^, turbulent ampli- 
tudes drop rapidly so that 11 decreases, resulting in the 
local maxima seen in Fig. [2J This suppression of 11 is due 
to linear stabilization as 7^ approaches 7^ from below. 
The momentum flux then increases monotonically when 

(c) 

lE > 1e ^^^ ^^ subcritical turbulence. The maxima in 
n may lead to a bifurcation in flow shear, discussed in 
the Conclusions. 

Prandtl number. The turbulent Prandtl number can 
be calculated from the values of Qi and 11. It is defined 
as Pr = Vi/xij where the turbulent thermal diffusivity, 
Xi, is given by Qi — ~XidTi/dr. Fig. |4] shows that Pr 
is approximately independent of k and only has strong 
dependence on •je for small values of 7^; [23]. This is 
despite the fact that both Vi and Xi individually have 
strong dependencies on k and je- For je ^ 0-4 the 
Prandtl number is close to unity, in good agreement with 
experimental measurements at low Mach numbers ^^. 

Stiff transport. A serious impediment to confinement 
is the sensitive dependence of the heat flux on small 
changes in k. This stiffness of the transport makes it 
difficult to increase k, and therefore the core tempera- 
ture, beyond the critical value, Kc, at which turbulence 
is excited. Recent experimental results indicate that stiff- 
ness, dQi/dn, may be reduced at low magnetic shear, s, 
and large values of 7^ [TT]. For the Cyclone base case 
considered here (s — 0.8), we find a complicated depen- 
dence of stiffness on flow shear, as illustrated in Fig. [5] 
At low values of 7^ (< 0.3), the critical k shifts to higher 
values, but the stiffness increases. This makes sense in- 
tuitively: one expects that for k/je — ^ 00, the heat flux 
will tend to the curve corresponding to je = 0. For this 







FIG. 5; Turbulent heat flux vs Ro/Lt for various values of 
the flow shear. Though there is some minimal relaxation of 
the profile stiffness for 0.3 < 7b < 0.8, the predominant effect 
of the shear is to shift the critical temperature gradient. 



ent. However, a detailed transport analysis is necessary 
to definitively observe such a bifurcation for experimen- 
tally relevant conditions. Finally, we note that the E x B 
suppression and PVG drive terms differ by a factor pro- 
portional to {qRo/r){B^/B) [5]. Thus our results, which 
should be qualitatively robust, may undergo considerable 
quantitative variation with changes in magnetic configu- 
ration. 
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to happen, the value of dQi/dK at Qi = must increase 
to compensate for the increase in Kc- 

For 0.3 < 7b < 7]^ , both Kc and the profile stiffness 
decrease. This is a result of a change in the nature of 
the linear instability, which is now driven by the PVG. 
The relaxation of stiffness is modest, and it only occurs 

(c) 

for K near Kc- When je > 1e i ^c initially shifts up- 
wards with increasing 7^, and stiffness increases for k, 
near Kc- However, for even larger je, both the stiffness 
and Kc decrease. This is to be expected, as the subcriti- 
cal turbulence is driven by the velocity, not temperature, 
gradient and thus has a weaker dependence on k. 

Conclusions. We have shown that PVG-driven tur- 
bulence exists at large flow shears, but only if the initial 
fluctuation amplitudes are sufficiently large. This has 
a number of potentially important implications. First, 
it indicates that linear stability analysis in insufficient to 
determine critical gradients in rotating plasmas. Further- 
more, it shows that the system can undergo hysteresis: 
The fluxes associated with a given flow shear and temper- 
ature gradient pair depend on the path taken to obtain 
that pair. As an example, consider the pair 7^; = 1, 
Rq/Lt = 10. From Fig. [2) we see that if Rq/Lt = 10 is 
obtained at je < Je before increasing 7s to unity, there 
will be large heat and momentum fluxes. In this case, 
there is an optimal flow shear for conflnement (7^; w 1) 
and a maximum attainable temperature gradient for a 
given power input. However, if 7^; = 1 is obtained at 
Rq/Lt ^ 7.5 before increasing Rq/Lt to ten, then the 
fluxes will be zero because the initial amplitude is not 
large enough to excite the subcritical turbulence. Thus 
it may be beneflcial to increase flow shear on a time scale 
short compared to the energy conflnement time. 

The existence of local maxima of the toroidal angular 
momentum flux provides the potential for bifurcations in 
^E and corresponding bifurcations in temperature gradi- 
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